library(tidyverse)

g2r_hg18 <- 106204113

g2r_hg19 <- 105737776

# ref for impute2
snp_legend <- read_delim('CRC-I/ref/ALL.chr14.integrated_phase1_v3.20101123.snps_indels_svs.genotypes.nosing.legend.gz')

# maf from TCGA-COAD --------
maf_tibble <- 'CRC-I/data/TCGA/COAD_maf/' |>
  list.files(pattern = 'maf.gz$',
           full.names = TRUE,
           recursive = TRUE) |>
  read_delim(comment = '#')

g2r_impute <- maf_tibble |>
  filter(Chromosome == 'chr14',
         Start_Position > g2r_hg19 - 2e6 | Start_Position < g2r_hg19 + 2e6) |>
  transmute(id = str_glue('{Chromosome}_{Start_Position}'),
            rsid = dbSNP_RS,
            position = Start_Position,
            allele = Reference_Allele,
            alt = Allele,
            individual = str_sub(Tumor_Sample_Barcode, end = 12))

g2r_pivot_indiv <- g2r_impute |>
  mutate(genotype = 1) |>
  pivot_wider(names_from = 'individual',
              values_from = 'genotype',
              values_fn = max,
              values_fill = 0,
              id_cols = position)


g2r_pivot_indiv |>
  transmute(chr = 'chr14',
            start = position,
            end = position) |>
  write_tsv('CRC-I/results/TCGA_COAD_ighg1_hg19.bed', col_names = FALSE)

# covert in UCSC liftOver... -------
g2r_hg18 <- read_tsv('CRC-I/results/TCGA_COAD_ighg1_hg18.bed',
                     col_names = c('chr','start','end','old','x'))

g2r_hg18_to_imp <- g2r_hg18 |>
  transmute(hg19 = str_extract(old, '(?<=-).+') |> as.numeric(),
            hg18 = start) |>
  filter(hg18 %in% snp_legend$position) |>
  left_join(g2r_impute, join_by(hg19 == position)) |>
  mutate(genotype = 1) |>
  pivot_wider(names_from = 'individual',
              values_from = 'genotype',
              values_fn = max,
              values_fill = 0) |>
  filter(rsid != 'novel') |>
  relocate(id, rsid) |>
  select(-hg19) |>
  mutate(across(6:last_col(), \(x)case_when(x == 1 ~ '1,0,0',
                                            x == 0 ~ '0,0,0'))) |>
  separate_wider_delim(6:last_col(), delim = ',', names_sep = '_')

# gens and .strand for impute2 --------
g2r_hg18_to_imp |>
  write_delim(col_names = FALSE, 'CRC-I/results/TCGA_COAD_ighg1_to_impute.gens')

g2r_hg18_to_imp |>
  transmute(position = hg18, strand = '+') |>
  write_delim(col_names = FALSE, 'CRC-I/results/TCGA_COAD_ighg1_to_impute.strand')

# impute2 output ---------
res_impute <- read_delim('CRC-I/results/TCGA_D5_6540_ighg1.impute2',
                         col_names = FALSE)

patients <- colnames(g2r_hg18_to_imp) |>
  str_subset('TCGA.+_1$')

predict_genotype <- res_impute |>
  filter(X3 == g2r_position) |>
  select(X6:last_col()) |>
  as.numeric() |>
  matrix(ncol = 3, nrow = 29, byrow = TRUE) |>
  as_tibble() |>
  set_names(c('GG','GR','RR')) |>
  add_column(patient = patients) |>
  pivot_longer(1:3, names_to = 'genotype', values_to = 'probability')

predict_genotype |>
  ggplot(aes(patient, probability, fill = genotype)) +
  geom_col()
